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Abstract 

What is the relationship between the macroscopic pa- 
rameters of the constitutive equation for a granular 
soil and the microscopic forces between grains? In or- 
der to investigate this connection, we have simulated 
by molecular dynamics the oedometric compression 
of a granular soil (a dry and bad-graded sand) and 
computed the hypoplastic parameters h s (the gran- 
ular skeleton hardness) and rj (the exponent in the 
compression law) by following the same procedure 
than in experiments, that is by fitting the Bauer's 
law e/eo = exp(— (3p/h s ) n ), where p is the pressure 
and eo and e are the initial and present void ratios. 
The micro-mechanical simulation includes elastic and 
dissipative normal forces plus slip, rolling and static 
friction between grains. By this way we have explored 
how the macroscopic parameters change by modify- 
ing the grains stiffness, V ; the dissipation coefficient, 
7„; the static friction coefficient, /i s ; and the dynamic 
friction coefficient, fik ■ Cumulating all simulations, 
we obtained an unexpected result: the two macro- 
scopic parameters seems to be related by a power 
law, h s = 0.068(4)?7 9 - 88 ( 3 ). Moreover, the experimen- 
tal result for a Guamo sand with the same granulom- 
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etry fits perfectly into this power law. Is this relation 
real? What is the final ground of the Bauer's Law? 
We conclude by exploring some hypothesis. 

1 Introduction 

Granular media are present everywhere. Examples 
go from common salt at the kitchen to planetary ice 
rings. In tons amount, granular media are the sec- 
ond most used materials, after water [TJ [5]. One of 
the most interesting examples of granular media are 
soils. A good understanding of the behavior of soils 
under several conditions is determinant in terms of 
engineering design, building planning and construc- 
tion processes. Furthermore, soils and granular me- 
dia represents a new paradigm in physics, and simula- 
tions by computers have turned out to be an excellent 
tool to gain deep insight in their behavior. 

Traditionally, two main streams have been used to 
understand soils O IH [5] . On one hand, engineers 
propose macroscopic constitutive relations in order 
to reproduce the deformation (or deformation rate) in 
terms of the strain (or strain rate) of the soil. Many 
formulations, like viscoplasticity, plasticity and hy- 
poplasticity have been successful to reproduce the 
experimental results. For instance, the hypoplastic 



model has been very useful to reproduce the experi- 
mental behavior of dry sands under monotonic loads. 
This model uses eight parameters to characterize the 
soil, and all of them can be obtained from element 
tests. On the other hand, physicists try to under- 
stand the soil behavior as the global result of micro- 
scopic forces between grains El H IE] • This global 
behavior is usually investigated by means of com- 
puter simulations [9]. To determine the relationship 
between the macroscopic parameters of the constitu- 
tive equation for a granular soil and the microscopic 
forces between grains is one of the main questions in 
the field. 

In this work we explore this connection for the case 
of a low polydisperse (bad-graded) dry sand when 
modeled by the hypoplastic theory. In particular, we 
want to investigate the dependence of two hypoplas- 
tic parameters that are obtained from the oedomet- 
ric test on the soil: the granular skeleton hardness, 
h s , and the exponent in the compression law, 77, in 
terms of the parameters governing the microscopic 
interactions between grains. For this purpose we per- 
form three-dimensional discrete element simulations 
of this element test on a dry sand of spherical grains 
for several combinations of four microscopic parame- 
ters, namely: the stiffness of the grains V, the normal 
damping coefficient 7„, the kinetic friction coefficient 
/ifc and the static friction coefficient /j, s . The soft- 
ware, developed by us, also includes rolling forces and 
torques, and therefore it is able to reproduce global 
reorganizations by rolling. The microscopic parame- 
ters are varied around those for a Guamo sand [TU] ■ A 
simulated oedometric test is performed for each set of 
microscopic parameters, and the two hypoplastic pa- 
rameters are measured from the simulation for each 
case. Finally, all simulations are put together in order 
to intend (if any) an empirical relation between the 
two macroscopic parameters, and the results are com- 
pared with the experimental values of the hypoplastic 
parameters for a Guamo sand. Section [2] introduces 
the microscopical forces included. Section [3] shows 
the integration algorithms employed. The simulated 
oedometric test are performed and analyzed in Sec. 
[4j Finally, Section [5] summarizes the conclusions and 
introduces suggestions for further research. 



2 The microscopic forces and 
torques 

The experimental system we want to simulate is a 
dry sand with very low polydispersity (between 0.85 
and 1.15mm in diameter). We model the grains as 
spherical particles in three dimensions. The torques 
and forces between grains act in normal and tangen- 
tial directions and dissipate energy on both of them. 
In the following, the subindexes i and j represent 
particle i and j, respectively, and the subindex ij de- 
notes relative quantities, (see Figure [I] and table [I] 
for details). 
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Figure 1: Three-dimensional vector quantities for two 
spherical grains in contact. 

The total force in normal direction is given by 
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where the first term is the Hertz elastic force and 
the second one is a normal damping force that repro- 
duces the experimental effect of the restitution coeffi- 
cient [Hj [12] . The variable hij = Ri+Rj — \fi — fj | is 
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tween grains depend on the relative position and mo- 
tion: Furthermore, in order to compute the torque 
j\ye need to define the location of the contact point 
i.e. the point where the forces are applied. It is 
usual to define c to be in the middle of the contact 
surface, despite the polydispersity of the sample, but 
Young modulus of grain i one can show by simple geometrical arguments that 
Poisson modulus of grain i this point c is located at [15] 
Normal unitary vector 



Radius of particle 
Position of the contact point i 
Corrected Radius of particle i 

Stiffness of grain i 



Normal tangential vector 
Normal component of vector B 
Tangential component of vectorrfie value zujy 
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-hij/2 requires Ri=Rj. 



Reduced mass In order to compute the sliding friction force and 

Rolling angle torqu e, we check out the value of the tangential rela- 



Table 1: 
model 



Parameter definitions of the microscopic 



the inter-penetration distance between particles. The 
normal unitary vector is computed as hij — (fi — 
rj)/\?i — In all cases but static friction, the tan- 
gential unitary vector is computed as Uj = 
where v\j is the tangential component of the relative 
velocity In addition, in the normal direction we have 
a torque that slows down the relative angular rota- 
tions on the normal direction. As an alternative to 
the traditional Cundall-Strack approach O E] , we 
derived a new expression for this torque [15] . It reads 



This expression can be deduced in the following way: 
when two grains are touching each other and have a 
non-null relative angular velocity on the normal di- 
rection, the contact surface of one grain rotates rel- 
ative to the other one. For each area element on 
the contact surface we assign a kinetic friction force 
that opposes to the relative motion. The net sum of 
this forces is zero, since they cancels in pairs, but the 
torque does not cancel. By summing up all torque 
contributions on the contact surface, we obtain ([2]). 
This expression depends on the reduced radius and 
the penetration depth as a consequence of the geome- 
try of the contact surface without additional assump- 
tions. 

On the tangential direction, the forces acting be- 



tive sliding velocity at the contact point 
Vpij = v pi - v pj = Vi- vj - {R'.fii + R'jti 



(4) 



If this velocity is different from zero, the particles are 
sliding and we apply the following force and torque: 



Ft. 
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= -R'h,i x Ft 



(5) 
(6) 



If the sliding velocity is almost null (in fact, less than 
a small predefined velocity) we compute the objective 
rolling velocity [T3] as 
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where is the reduced effective radius. If this ve- 



locity is non-null, the particles are rolling. By ex- 
tending the rolling model to three dimensions and 
applying on two soft spherical bodies we obtain the 
following expression for the rolling friction force and 
torque [T5] . 
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where the rolling constraint is clear. Finally, if the 
particles are not sliding and not rolling over each 
other, they are in relative rest. But, if there is a tan- 
gential relative force, the particles could not stay in 
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rest unless a tangential static friction force is present. 
In this case the tangential unitary vector iy is com- 
puted from the tangential force as = F t /\F t \. 
For the static friction force we apply a simplified 
model [TBI. 
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With the model presented so far we were able to 
reproduce complex behaviors as sliding or rolling and 
dissipation on both normal and tangential directions. 
We tested each force implementation by making par- 
ticular simulations for each of the following cases: 
sliding, rolling, only static friction and normal dis- 
sipation. For each case the model worked properly 
and also for the case with all the interactions turned 



3 Integration methods 

The simulation method we used is Molecular Dynam- 
ics (MD), also known as Discrete Element Method 
(DEM). In MD, the time evolution is traced on dis- 
crete time steps of size St. The size of St depends on 
the mechanical and geometrical properties of the sys- 
tem and is typically of the order 0(St) ~ 10~ 5 — 10 
s. The MD computes the next positions and veloci- 
ties (or the next angular orientation and angular ve- 
locities) by solving the second Newton law for each 
particle in terms of the current positions, velocities 
and forces (or the current orientation, angular veloc- 
ities and torques). This needs not only a model for 
the forces and torques, as given in Section[2j but also 
an integration algorithm. 

There exists many different integration algorithms 
for the translation and rotation variables. The choice 
of one algorithm over another depends on the forces, 
the system size and the total simulation time. In 
order to choose one integration algorithm for the 
translational variables we investigated the conserva- 
tion of energy in a system of 50 particles colliding 
with the Hertz force in a rectangular box. The im- 
plemented translational integration algorithms were: 
Verlet [13 US], Leap-Frog [HI HI], optimized Velocity 
Verlet [HI [20] and a fifth-order predictor-corrector 



method [T7]- The optimized velocity Verlet method 
is written as 

Ri = R(t) + V(t)(,St, (11) 

Vl = V(t) + —F[R!]St/2, (12) 

m 

4 = ^1+^(1-20^, (13) 

V(t + St) =Vi+ -F[R 2 ]St/2, (14) 
m 

R(t + St) =R 2 + V(t + 5t)£5t, (15) 

where R, V and F represents the particle position, 
velocity and force, respectively. The parameter 
£ ~ 0.193183325037836 minimizes the total trunca- 
tion error of the algorithm |19j . For each time step, 
two computations of the forces are needed, but even 
doubling the time step makes the errors three times 
smaller than in the original Velocity Verlet algorithm 
(£ = 0) . Figure [2] shows the variation of the total 
mechanic energy (SE 2 ) 1 / 2 = ((E 2 ) — (E) 2 ) 1 / 2 , aver- 
aged over the total simulation time, as a function of 
the time step St. It is clear that for large St the best 
algorithm is the optimized velocity Verlet method. 
This is the algorithm we chose for the translational 
motion. 
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function of the time step St for different translation 
integration algorithms. 

The simulation of spatial rotations is more com- 
plicated. It is well known that the Euler angles 
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can keep the information of a rigid body on space, 
but there are some numerical instabilities on using 
them [TTl [18] . For this reason, the Euler angles are 
replaced by unitary quaternions, denoted here by q, 
and the numerical problems are solved. The main 
disadvantages on using quaternions are the absence 
of high order integration algorithms and the need to 
normalize the quaternion at each time step. These 
problems are solved in a new formulation proposed 
by Omelyan [3T] of the Leap- Frog method for quater- 
nions that preserves the quaternion unit norm, de- 
spite the size of the discrete time step. 

In the original paper, the Omelyan algorithm is 
formulated in the y convention of the Euler angles. 
We rewrote it on the x convention, 



Q(uj)q 



(16) 

while the new angular velocities and quaternions are 
computed as 



initial vertical pressure po is applied on the top wall 
and, after some time (typically 10 minutes), the sam- 
ple relaxes to a new void ratio. Then, the pressure 
is doubled, the sample relaxes and a new void ratio 
is measured. The procedure repeats until reaching 
some maximal pressure before the particles crush. 

The curve so obtained (figure [3]) can be fit by the 
Bauer empirical equation |22] 
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where the superscript b represents the body fixed ref- 
erence axes, I is the identity matrix on Br, Q(uj) is 



the matrix defined on ( 16 ) and the angular velocity 
tu is computed in t + St/2. 

4 Simulation of the oedometric 
test 

Several element tests are performed on a soil in order 
to measure its macroscopic parameters. An oedomet- 
ric test [7] consists on filling a metallic cylinder with 
a sample of the soil and measuring the initial void ra- 
tio eo, that is the ratio between the volume occupied 
by voids over the total volume of the sample. All 
cylinder walls but the top one are fixed. Then, an 



e 
eo 



exp 



(19) 



by using, for instance, the Levenberg-Marquardt al- 
gorithm [23] . This gives the two hypoplastic param- 
eters h s and rj. Figure [3] shows the oedometric curve 
for a Guamo sand with grains between 0.85mm and 
1.15mm diameter; we obtained ft, s [MPa] = 19.1(1.1) 
and n = 0.57(1). 
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Figure 3: Experimental oedometric test results on a 
dry Guamo sand in (top) linear and (bottom) log- 
scale on the applied vertical stress. 
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The same experimental procedure has been im- 
plemented in the simulations. For each set of mi- 
croscopic parameters, the grains are randomly dis- 
tributed inside the cylindrical container; they fall 
down by gravity and collide each other and with the 
walls until rest. Then, the top wall falls down with 
the initial pressure po and the oedometric test starts, 
following the same procedure than in the experimen- 
tal test. We try the simulation to be as similar to the 
experiment as possible: The dimensions of the con- 
tainer, the polydispersity of the sample and others 
parameters have almost the same experimental val- 
ues. The only difference is in the number of grains: 
around 5000 for the experiment and 280 for the sim- 
ulation (because of hardware limitations) . At a first 
glance it appears as a very drastic reduction, but the 
results were very closed to the experimental data, and 
it gives us some confidence on the procedure. More- 
over, current simulations with more than 1000 par- 
ticles shows similar behaviour to the present ones. 
The typical computational time for the simulation of 
a complete oedometric test was about 52 hours on 
a Pentium IV 3.0 GH machine with 4GB of RAM. 
The compiler used was gnu/g++. That was done for 
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Table 2: Macroscopic parameters h s and r\ for differ- 
ent values of the microscopic parameter V (top) and 
Hk (bottom). 

Two typical curves for the simulation of the oedo- 



1 Currently, by using the same theoretical model and some 
special compiler flags, we are able to simulate the oedometric 
test in about three days with 1200 particles. 
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Figure 4: Void ratio versus vertical stress for various 
values of the microscopic static friction coefficient /i s . 



metric compression are shown in Figure [4j Similar 
curves are obtained for the other parameters by fol- 
lowing the same procedure as before (see Tables [2] 
and [3] for details). 

We noted that, in all cases, the two hypoplastic 
parameters behave in a non-independent way. When 
one parameter increases, the other one decreases. 
This suggest a hidden relationship between them. 
In order to investigate this relation, we plot h s as 
a function of 77 for all sets of microscopic parame- 
ters, as shown in figure [5j It reveals and unsuspected 
power law relationship between the two hypoplastic 
parameters: h s —0.068(A)ri~ 9 88 ^. Moreover, the ex- 
perimental values lay on the line. This kind of rela- 
tion could suggest a possible reduction in the number 
of macroscopic parameters of the hypoplastic theory, 
but the experimental values for different sands with 
different granulometries do not lay on a single power 



6 



7 „[s- 1 m- 1 /2] 


*7 




/i s [MPa] 


5.0el 


287(1) xlO 


-3 


250(4) xlO 2 


5.0e2 


29(4) xlQ- 


-2 


16(8) xlO 3 


5.0e3 


32(2) xlQ" 


-2 


8(1) xlO 3 


5.0e4 


27(4) xlO" 


-2 


2(1) xlO 4 


5.0e6 


26(3) xlO" 


2 


4(2) xlO 4 


5.0e7 


20(4) xlO" 


-2 


7(6) xlO 4 



Us 



h s [MP£ 



0.1 
0.2 
0.4 
0.6 
0.7 
0.8 



26(2) xlO 
30(2) 
24(1) 
20(2) 
19(3) 
23(1) 



xlO- 2 
xlQ- 2 
xlQ- 2 
xlQ- 2 
xlQ- 2 



32(7) xlO 3 
7(2) xlO 3 

4(1) xlO 4 
5(3) xlO 5 

10(8) xlO 5 

5(2) xlO 5 



Table 3: Macroscopic parameters h s and rj for differ- 
ent values of the microscopic parameter j n (top) and 
H s (bottom). 
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Figure 5: h s as a function of r\ for different values of 
the microscopic parameters 



law. Let us point out that we have varied all mi- 
croscopic parameters (V, 7„, /i& and /i s ), but the 
granulometry. It is well known [24 that r\ strongly 
depends on the granulometry, but once this variable 
is fixed, a power law could appear. This result should 
be validated by future experiments and simulations. 



5 Conclusions 

Hereby we have simulated the oedometric compres- 
sion of a bad-graded Guamo sand in order to esti- 
mate the hypoplastic parameters h s and 77, and we 
have compared with experimental results. In partic- 
ular, we have investigated the changes on these two 
macroscopic parameters by varying four microscopic 
parameters: the stiffness V, the normal damping con- 
stant 7„, the static friction coefficient fi s and dynamic 
friction coefficient /_tfc. By doing so, we found an unex- 
pected power-law relationship between h s and 77: for 
our case, h s cx jy- 9 - 88 ( 3 ). Moreover, the experimental 
values for the Guamo sand with the same granulome- 
try (random-distributed diameters between 0.85 and 
1.15mm) lies on the same curve. This suggests us 
that these two parameters may be replaced by other 
two, more related with the microscopic world: one 
reflecting the granulometry and another one reflect- 
ing the strength of the microscopic interactions. This 
power law, of course, must be confirmed by a broader 
set of experiments, but the possibility is very promis- 
ng. The experimental confirmation of these kind of 
■elationships and the possible definition of these two 
rew parameters are subjects of present research. 

The simulation performed is 3D and the soft- 
ware developed by us includes some of the state- 
the-art algorithms for spatial translation and ro- 
ations. Moreover, the microscopic force model in- 
dudes rolling among the more traditional elastic and 
lamping normal interactions and sliding and static 
rictional forces. The rolling process allows global 
ttjssipation and long-range reorientations without re- 
quiring a big spatial reconfiguration of the sample, 
that is without changing the macroscopic void ratio, 
and its global effect on the macroscopic parameters 
is very interesting to explore. For example, it would 
be possible to perform simulations with and without 
the rolling force in order to get insight into its actual 
role on the soil properties. This will be a topic of 
future work. 

Micro-mechanical simulations are powerful tools 
for the investigation on the microscopic origin of 
macroscopic behavior of soils. Furthermore, these 
computer experiments have shown to be able to ap- 
proximate the experimental parameters for the soil. 
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We hope this work will raise new questions regard- 
ing the microscopic origin of the constitutive param- 
eters and gives starting points for a redefinition of 
the macroscopic parameters, more related with the 
microscopic world. 
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